// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_HOUSEHOLDER_H
#define EIGEN_HOUSEHOLDER_H

namespace Eigen {

namespace internal {
    template <int n> struct decrement_size
    {
        enum
        {
            ret = n == Dynamic ? n : n - 1
        };
    };
}  // namespace internal

/** Computes the elementary reflector H such that:
  * \f$ H *this = [ beta 0 ... 0]^T \f$
  * where the transformation H is:
  * \f$ H = I - tau v v^*\f$
  * and the vector v is:
  * \f$ v^T = [1 essential^T] \f$
  *
  * The essential part of the vector \c v is stored in *this.
  * 
  * On output:
  * \param tau the scaling factor of the Householder transformation
  * \param beta the result of H * \c *this
  *
  * \sa MatrixBase::makeHouseholder(), MatrixBase::applyHouseholderOnTheLeft(),
  *     MatrixBase::applyHouseholderOnTheRight()
  */
template <typename Derived> EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
{
    VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size() - 1);
    makeHouseholder(essentialPart, tau, beta);
}

/** Computes the elementary reflector H such that:
  * \f$ H *this = [ beta 0 ... 0]^T \f$
  * where the transformation H is:
  * \f$ H = I - tau v v^*\f$
  * and the vector v is:
  * \f$ v^T = [1 essential^T] \f$
  *
  * On output:
  * \param essential the essential part of the vector \c v
  * \param tau the scaling factor of the Householder transformation
  * \param beta the result of H * \c *this
  *
  * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(),
  *     MatrixBase::applyHouseholderOnTheRight()
  */
template <typename Derived>
template <typename EssentialPart>
EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholder(EssentialPart& essential, Scalar& tau, RealScalar& beta) const
{
    using numext::conj;
    using std::sqrt;

    EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
    VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size() - 1);

    RealScalar tailSqNorm = size() == 1 ? RealScalar(0) : tail.squaredNorm();
    Scalar c0 = coeff(0);
    const RealScalar tol = (std::numeric_limits<RealScalar>::min)();

    if (tailSqNorm <= tol && numext::abs2(numext::imag(c0)) <= tol)
    {
        tau = RealScalar(0);
        beta = numext::real(c0);
        essential.setZero();
    }
    else
    {
        beta = sqrt(numext::abs2(c0) + tailSqNorm);
        if (numext::real(c0) >= RealScalar(0))
            beta = -beta;
        essential = tail / (c0 - beta);
        tau = conj((beta - c0) / beta);
    }
}

/** Apply the elementary reflector H given by
  * \f$ H = I - tau v v^*\f$
  * with
  * \f$ v^T = [1 essential^T] \f$
  * from the left to a vector or matrix.
  *
  * On input:
  * \param essential the essential part of the vector \c v
  * \param tau the scaling factor of the Householder transformation
  * \param workspace a pointer to working space with at least
  *                  this->cols() entries
  *
  * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), 
  *     MatrixBase::applyHouseholderOnTheRight()
  */
template <typename Derived>
template <typename EssentialPart>
EIGEN_DEVICE_FUNC void MatrixBase<Derived>::applyHouseholderOnTheLeft(const EssentialPart& essential, const Scalar& tau, Scalar* workspace)
{
    if (rows() == 1)
    {
        *this *= Scalar(1) - tau;
    }
    else if (tau != Scalar(0))
    {
        Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols());
        Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows() - 1, cols());
        tmp.noalias() = essential.adjoint() * bottom;
        tmp += this->row(0);
        this->row(0) -= tau * tmp;
        bottom.noalias() -= tau * essential * tmp;
    }
}

/** Apply the elementary reflector H given by
  * \f$ H = I - tau v v^*\f$
  * with
  * \f$ v^T = [1 essential^T] \f$
  * from the right to a vector or matrix.
  *
  * On input:
  * \param essential the essential part of the vector \c v
  * \param tau the scaling factor of the Householder transformation
  * \param workspace a pointer to working space with at least
  *                  this->rows() entries
  *
  * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), 
  *     MatrixBase::applyHouseholderOnTheLeft()
  */
template <typename Derived>
template <typename EssentialPart>
EIGEN_DEVICE_FUNC void MatrixBase<Derived>::applyHouseholderOnTheRight(const EssentialPart& essential, const Scalar& tau, Scalar* workspace)
{
    if (cols() == 1)
    {
        *this *= Scalar(1) - tau;
    }
    else if (tau != Scalar(0))
    {
        Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace, rows());
        Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols() - 1);
        tmp.noalias() = right * essential;
        tmp += this->col(0);
        this->col(0) -= tau * tmp;
        right.noalias() -= tau * tmp * essential.adjoint();
    }
}

}  // end namespace Eigen

#endif  // EIGEN_HOUSEHOLDER_H
